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Lagrange-mesh calculations in momentum space 
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The Lagrange-mesh method is a powerful method to solve eigenequations written in configuration 
space. It is very easy to implement and very accurate. Using a Gauss quadrature rule, the method 
requires only the evaluation of the potential at some mesh points. The eigenfunctions are expanded 
in terms of regularized Lagrange functions which vanish at all mesh points except one. It is shown 
that this method can be adapted to solve eigenequations written in momentum space, keeping the 
convenience and the accuracy of the original technique. In particular, the kinetic operator is a 
diagonal matrix. Observables in both configuration space and momentum space can also be easily 
computed with a good accuracy using only eigenfunctions computed in the momentum space. The 
method is tested with Gaussian and Yukawa potentials, requiring respectively a small or a great 
mesh to reach convergence. 
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I. INTRODUCTION 

There are few three-dimensional problems in quantum mechanics which allow a complete analytical solution for any 
value of the orbital angular momentum (the S'-wave channel is very similar to a simpler one-dimensional equation) 
[l|. So, numerous methods have been developed to solve numerically with a high accuracy the eigenvalue equations 
associated with various systems. Among these techniques, the Lagrange- mesh method (LMM), which is especially 
easy to implement, can produce very accurate results. First created to compute eigenvalues and eigenfunctions of 
a two-body Schrodinger equation [2-01' it has been extended to treat semirelativistic Hamiltonian |8l-[lll|. The trial 
eigenstates are developed in a basis of particular functions, the Lagrange functions, which vanish at all mesh points, 
except one. Once the potential is known in the configuration space, its matrix elements are simply the potential values 
at the mesh points, if they are computed with an associated Gauss quadrature. At first sight, this method could look 
like a discrete variational method, but this is absolutely not the case since the eigenfunctions can be computed at 
any position. Because of the use of the Gauss quadrature scheme, the method is not variational but the results are 
too a large extent independent of the sole nonlinear parameter of the method fixing the physical scale of the system. 
Generally, a great accuracy can be reached with a small mesh jl2 |. 

At the beginning, the LMM was developed in the configuration space. Recently, it has been shown that the Fourier 
transform of the eigenfunctions computed in the configuration space can easily be obtained with a good accuracy in the 
physical domain of the momentum space [1.3-|. Moreover, observables in this space can easily be computed with a good 
accuracy using only matrix elements and eigenfunctions in the configuration space. But, for some particular problems, 
it can be preferable to work in the momentum space. This is the case when the potential presents discontinuities in the 
configuration space [14] or when the potential is given in the momentum space. In this last case, if it is possible to use 
the LMM by computing first the Fourier transform of the potential, we will show here that the LMM can be adapted 
to solve the eigenequations directly in momentum space. Observables in both configuration and momentum spaces 
can also be computed with a good accuracy using only eigenfunctions computed in momentum space. Moreover, 
we will show that the new LMM can provide these types of data very efficiently and very easily, using again the 
fundamental properties of the Lagrange functions. Let us note that the method presented here relies on a mesh of 
points built with the zeros of a Laguerre polynomial, but a general procedure for deriving other Lagrange meshes 
related to orthogonal or non-orthogonal bases has also been developed [la] . 

Several methods exist to solve eigenequations in momentum space. For instance, iterative procedures have been 
developed [16lll7| . They are quite accurate but resort finally to numerical integrations on a mesh. Direct computations 
on a mesh are easier to carry out, but they require very large mesh if a good quadrature rule is not used |14| . As we 
will see, the LMM in momentum space is very easy to implement and can also give accurate results. In order to fix the 
notations, the eigenequations in momentum space are presented in Sec. |TT1 The LMM adapted in momentum space 
is described in Sec. Illll with some details in order that the paper is self-contained. Test calculations are presented in 
Sec. II VI for two potentials with very different properties of convergence. Some concluding remarks are given in SecFVl 

II. EIGENEQUATIONS IN POSITION AND MOMENTUM SPACES 

Let us consider the following eigenequation (h = c = 1) 

[Tip') + Vir)]\^)^E\<l^), (1) 

in which the kinetic part and the potential depend respectively on the relative square momentum p^ and on the radial 
distance r = \r\ between the particles. The wavefunctions in configuration space 4)r{r) = {'r\4') and in momentum 
space (jipij)) = {p\4') can be written using the spherical representation 

<^,(r)=7^„^(r)l^,„(f), (2) 

C^j,{p)=Vnl{p)Yi.,n{P). (3) 



where Yim{x) = i^Yim{x) is a modified spherical harmonic [18|, x = x/\ 
distance. These functions are linked by the following Fourier transforms [l9| 



and n is the number of nodes at finite 



Mp) = J^^ J Mr) e-^^^-Ur, (4) 

^r (r) = j^^ J 4>p (p) e+^P-~ dp. (5) 



These equations lead to [ic 



Vni{p) = i-iy^^J nniir)jiipr)r^dr, (6) 



7^„^(r) = (-l)^/- / Vni{p)ji{pr)pUp, (7) 



where ji (x) is a spherical Bessel function J20| . 

Written in the momentum space, ([IJ takes the following form 



T{p^) Mp) + J Vft{p~ p') Mp') dp' = E (Ppip) (8) 

with Vft(p — p')i the Fourier transform of V{r), given by 

VMp- p') = J^ J V{r) e-^(P'P'y^^dr. (9) 

This potential is a continuous function of the momentum, even if parts of the interaction in configuration space present 
discontinuities. One can think of square well or Dirac delta function (repulsive only in a 3D-space). As the potential 
depends only on r, we have Vft{p~p') — ^ft(|p — p'|) and ^ becomes [2l| 

1 f°° 
VFT{k) = —-TY V{r) shi{kr)rdr. (10) 
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Using the standard decomposition of a radial function [18] , the eigenvalue equation (|5]) takes the form of an integral 
equation for the wavefunction Vni (p) 



Tip^) Vnlip) + I Vlip,p') Vnl{p')p'^ dp' = EVnlip) (H) 







with the partial potentials 



Vi{p,p')=2'K I Pilt) Vft (\/p^+p'^ - 2pp't\ dt. (12) 

The Legendre polynomial Pi{t) depends on the variable t = p ■ p' . 
For a Schrodinger equation, the kinetic operator is given by 

T{p )^— with /^ = ■ 13) 

2/i nil + "12 

and the eigenvalue E is the binding energy of the system. In a spinless Salpeter Hamiltonian, the kinetic operator 
takes the following form 



T{p') = ^p^ + ml + ^p^+ml (14) 

and the eigenvalue E is the mass of the system. This kind of Hamiltonian is sometimes denoted semirelativistic since 
it is not a covariant formulation. This equation can be considered as a Schrodinger equation with its nonrelativistic 
kinetic part replaced by a relativistic counterpart. More rigorously, it is obtained from the covariant Bethe-Salpeter 
equation [22] with the following approximations: Elimination of any dependences on timelike variables and neglect of 
particle spin degrees of freedom as well as negative energy solutions 1231 . T he spinless Salpeter Hamiltonian is often 
used in hadronic physics to study bound states of quarks or gluons |24l - [27f . 

Within this formulation, the action of the kinetic operator is just an ordinary multiplication. So, nonrelativistic 
and semirelativistic systems are computed with the same manner. Moreover, more complicated kinetic parts, with 
momentum-dependent masses |28l430l | , can be equally treated. Though the formulations in the configuration and mo- 
mentum spaces are completely equivalent, this does not mean that the technical difficulties to solve the eigenequations 
are the same in both spaces. 

If the potential is known in the configuration space, ([TU| and ([T^ allows the computations of the partial potentials 
Vi(j),p'). VpT{k) can also be directly obtained from a physical theory, like a field theory naturally written in the 
momentum space. One can think of effective potentials obtained from Feynman diagrams, for instance. It is not 
always possible to obtain an analytical form for (I12p . but such a numerical integration can be rapidly and accurately 
performed. 



III. METHOD IN MOMENTUM SPACE 

A. Lagrange functions 

The LMM relies on the existence of a iV-point mesh {xi}, which is associated with an orthonormal set of N 
indefinitely derivable functions fj{x), called the Lagrange function 2-4]. Each function fj{x) satisfies the Lagrange 
conditions, 

f^{x,)^x;'^'s,,, (15) 

that is to say it vanishes at all mesh points except one. The Xi and Xi are respectively the abscissae and the weights 
of a Gauss quadrature formula 



POO 



N 

g{x)dx ^'^Xkgixk). (16) 

fe=i 



Several quadratures are possible, but, as we work with wavefunctions depending only on the module of a variable, 
we consider the case of the Gauss-Laguerre quadrature whose domain of interest is [0, oo[. The Gauss formula ([T5)) is 
exact when g{x) is a polynomial of degree 2N — 1 at most, multiplied by exp(— a;). The mesh points Xi are the zeros 
of a Laguerre polynomial of degree N: L^(xi) = |2|- These zeros can be determined with a high precision with 
usual methods to find the roots of a polynomial [3l| (the Mathematica expression Root does efficiently the job) or as 
the eigenvalues of a particular tridiagonal matrix [32] . The weights can be computed by the following formula [12j | 

N 

InA, =a;j -lna;j + 21nr(7V + l) - ^ ln{x,~Xjf. (17) 

Wavefunctions pVniip) vanish at the origin. As it is not the case for the original Lagrange functions, it is preferable 
to use the regularized Lagrange functions whose explicit form is given by fi{x) = {x/xi)fi{x), that is to say [3|, |6| 

Mx) = i^iyx;'^^x{x - x,y^LN{x) eM-^f^)- (is) 

Such a function satisfies ((TS]) and /i(0) = (see Eq. [5D]). 

B. Matrix equation 



The use of the LMM in configuration space is described in [2t8|, |13| . We will present here the formulation for the 
integral equation (|lll) within the LMM. The idea is to expand the wavefunction Vniiv) with the regularized Lagrange 
functions in such a way that a trial state \il)) is written 

I^)=EQI/.) with (jJ\f^)^lM^YUp)- (19) 

This formula is identical to formula (6) in [ij] with just the replacement of the variable r by p. The coefiicients Cj 
are linear variational parameters and the scale factor /i is a nonlinear parameter aimed at adjusting the mesh to the 
domain of physical interest. The parameter h has here the dimension of a momentum {h plays the same role in [l3| 
but has the dimension of a distance). Contrary to some other mesh methods, the wavefunction is also defined between 
mesh points by (|18p and (fTO)) . For a good value of h and a sufficiently high value of TV (see sec. lIVp . the function 

vM = J:c,^-^ (20) 

can be a good approximation of the exact function Vni{p)- Note that an advantage of the LMM is that the value for 
h has not to be known with a great accuracy. It is sufficient that this value be located within a given interval, as we 
will check in the following. 



Basis states \fi) built with the regularized Lagrange functions are not exactly orthogonal. But, at the Gauss 
approximation, we have {fj\fi) = Sji. So, in the following, all mean values will be performed using the Gauss 
quadrature formula (IT6l) . Inserting expansion (|20| in ((TT|) gives 

N , . N ^ f ( \ 

T{h^x^) Y, C, ^^^ + J2 ^1 h"" V^o^o Vi{hx,,hx) = eY,Cj ^^^, (21) 

where x = p/h is a dimensionless variable. We can now multiply this equation by x fj{x) and integrate on [0, (X)[ with 
again the Gauss quadrature formula p6l) . Finally, we obtain 

N 

Y^ C-j \T{h^x^) S,j + h^ y/X~\]x, Xj Vi{hx,,hxj) - E 5.,j = 0. (22) 

The Hamiltonian matrix is symmetric since Vi{p,p') = Vi{p' ,p). A similar expression is obtained for calculations with 
the LMM in the configuration space for a nonlocal potential [5]. With the LMM in momentum space, the solution 
of a quantum equation reduces (as it is often the case) to the determination of eigensolutions of a given matrix. So, 
once the partial potentials are known, (P^ shows that this method is very easy to implement. The computations of 
the Fourier transform of the wavefunction and of observables is no more complicated, as presented in the following. 

In Sec. IIV[ we will study the convergence of the method as a function of the scale parameter h and the number 
of mesh points N for nonrelativistic an semirelativistic kinematics. Let us note that an automatic determination of 
h has been developed for the LLM in the configuration space [Ij, |3j|- But, the generalization of such a technique 
to the LLM in the momentum space is very difficult due to the nonlocal nature of the interaction in (fTTj) . We will 
cross-check our results by comparing eigenvalues and mean values of observables computed with the LMM in both 
configuration and momentum spaces. Moreover, a wavefunction computed in the momentum space will be compared 
with the Fourier transform of the corresponding wavefunction computed in the configuration space, and vice versa. 

C. Fourier transform 

It has been shown in 13] that a good approximation of a function Vni (p) can be obtained from the Fourier transform 
of the corresponding solution computed in the configuration space by the LLM. This can be performed by taking 
benefit of the very special properties of the regularized Lagrange functions. Similarly, a good approximation of the 
function TZni (r) can be obtained from the Fourier transform of a solution computed in the momentum space by the 
LLM. Using the Gauss quadrature rule P^ with the Lagrange condition ([T5|) . the integral ([7]) for a given trial state 
(pn| simply reduces to 

n ^ ^ 

^nz(r) = (-l)y-/i3/2^QVAlx,jK/ia;.r). (23) 

This formula is identical to formula (22) in [l^] with just the replacement of the variable r by p. This results from 
the similarity of dU and (O, and the choice of expansion (fT9)) . For a sufficiently high value of N (see sec. lIVp . TZniir) 
can be a very good approximation of the genuine function TZni{r) in the configuration space. Above a critical value 
of r, generally in the asymptotic tail of TZniir), this function can present large unphysical rapid oscillations. These 
oscillations do not develop in Vni{p), because they are damped by the rapid decrease of the regularized Lagrange 
functions. 

D. Mean values of momentum-dependent observables 

The mean value of the operator U{p) for a trial state \tp) is given by 

N 

imipM = E ^^ c-, {m{p)\f,). (24) 

Using the Lagrange condition ([Tsj) and the Gauss quadrature (fT6)) . this integral reduces to 

N 

{^\U{pm=Yc]U{hx,). (25) 



If U is the identity, we recover the normahzation condition as expected. As we will see in sec. lIVl a very good accuracy 
can be reached for the mean values {U{p)). 

E. Mean values of radial observables 

The mean value of the operator K{r) for a trial state \ip) is given by 

N 

{^\K{rM - J2 ^^ ^1 (MK{r)\fj). (26) 

The method to compute matrix elements {fi\K{r)\fj) relies on the fact that r^ = —V?. in the momentum space j34| . 
Let us first look at the matrix P whose elements are Pij — {fi\r^\fj)- With p^ . these matrix elements are given by 






^^. = 7:2 *^. + ^-::r^% ' (27) 



where 



Uo = j^ M^) (-^) /. W dx « -Xy'f^'ix,). (28) 

This expression is exact for some Lagrange meshes, but this is not the case for the regularized Laguerre mesh. An 
exact expression can easily be obtained (see the Appendix in Q). However, as shown in [4], it is preferable to use the 
approximation (l?fl) - (P5|) . The quantities tij are then easy to obtain and read [^ 

, _ I (— ) [XiXj) [Xi + Xj){Xi — Xj) {if^Jji (e)Q.\ 

^*^"\ (12x2)-i[4 + (4iV + 2)x,-a:f] [i^j). ^^^> 

If P^ is the diagonal matrix formed by the eigenvalues of P, we have 

P = SP^S-\ (30) 

where S is the transformation matrix composed of the normalized eigenvectors. Let us call K^ the diagonal matrix 
obtained by taking the function K{y/x) of all diagonal elements of P^ (remember that P is linked to the matrix 
elements of r^, not r). The numbers {fi\K{r)\fj) are well approximated by the elements of the matrix K obtained 
by using the transformation ^U\i: K = S K^ S^^. As we will see in sec. llVi a very good accuracy can be reached for 
the mean values {K{r)). 

IV. NUMERICAL RESULTS 
A. Gaussian potential 

We first consider a Gaussian potential whose Fourier transform is well known |2ll| 



V{r) = -a exp (-6^H) ^ VMk) = -^^^ST^e^p ^-— j . (31) 

Let us note that the asymptotic behavior of the potential is similar in the configuration and momentum spaces. Using 
(112), we can compute 



, „ a ( p2+p'2\[(Z5 ,fl\f2l-2k 



E2k-i[-^) + {-iy ''E2k-i(^ 



(32) 



where E,n{z) is an exponential integral [20] and [y\ means the integer part of y. This type of potentials is the prototype 
for a short-range interaction and can be used in many field of physics. 



With a nonrelativistic kinematics, we can write 



62 



where e is the solution of the dimensionless Hamiltonian 



H 



q' 



ge 



with g 



2 fia 
~62~' 



(33) 



(34) 



We choose to study the performance of the LMM in momentum space with the value g = 15. In this case, two bound 
states exist: {n,l) — (0,0) and (0,1). 

Let us look at the results obtained for H (p4|) by solving ((22)) . The variation of the eigenvalue e as a function of the 
nonlinear parameter h is presented for the two states in Fig. [TJ For a value of N as small as 20, a plateau is present 
with abrupt variations of the eigenvalue at the borders. The length of the plateau increases with the number N of 
mesh points. For the excited states, the convergence is usually slower than for the ground states: The plateau is less 
extended and the variations around the plateau are larger. Even with N = 10, very good results can be obtained 
provided h is taken in the small corresponding quasi-plateau (the fluctuation in the rapid variation). Note that the 
flatness of the plateau is a criterion obviously depending on the accuracy one wants to reach. It is shown on Fig.[5]that 
eigenvalues, corresponding to different values of h taken in plateaus, converge towards the same value as N increases. 
This convergence can be achieved from above or from below. All these results clearly show the non variational nature 
of the LMM. 
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FIG. 1. Eigenvalues e of (|22|) for the dimensionless Hamiltonian (|34|l with g = 15, as a function of h for three values of A'^. 
Left: Ground state. Right: Excited state with {n,l) — (0, 1). The scale for ordinates is the same for both graphics. 
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FIG. 2. Ground state eigenvalues e of (I22II for the dimensionless Hamiltonian p4|) with g — 15, as a function of A'^ for two 
values of h. For A'^ = 10 and h = 1.0, the value of e is —5.37859 far below the range of the graph. 



It is remarkable that, even for small values of N, the value of h must not be determined with a high precision. 
This is a great advantage for the LMM. Nevertheless, if h is too small, a significant part of the wavefunction is not 



8 

covered by the points of the mesh. When h is too large, all points of the mesh are located in the asymptotic tail 
of the wavefunctions and it is then impossible to obtain good eigenvalues. Let us note that the same behaviors are 
observed for the LMM in configuration space [2, Ij, Sl • 

Eigenvalues and some observables have been calculated with the LMM in both configuration and momentum spaces. 
Results are presented and compared in TableUfor the ground state only, but results are similar for the excited state. A 
very good accuracy can be obtained for the computation in momentum space, even for small values of N. For most of 
the physical problems, the required precision can probably be reached with N = 20. It seems that position-dependent 
observables converge more slowly than momentum-dependent observables. We have checked that the situation is the 
opposite for LMM calculations in configuration space. 

In principle, we must have (H) = e. This is not the case since momentum-dependent and radial observables are 
not computed with the same method (see Sec. IIII DJ and HITE)) . We can see in Table |T] that the difference between the 
two quantities is around the accuracy of the eigenvalues e. 



TABLE L Eigenvalues e and mean values of some observables for the ground state oi H = q^ + U{x) with U{x) = —ge ^ p4|) 
and g = 15, as a function of A*' for h = 0.5. Computations in momentum (Mom.) space are compared with accurate results 
obtained in configuration (Conf.) space for TV = 100 and h = 0.4. In order to check the mean values, {H) (which must be equal 
to e) is computed by (q^) -)- {U{x)). 
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-5.3776125307238 
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{^-') 


3.74063887622353 
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3.74063887622358 


(q') 


26.50642515647 


26.50643641212 


26.50642516641 


26.50642515646 


(x) 


0.7134620 
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0.7134650 


0.7134620 


{U{x)) 


-9.1182387832920 


-9.1182424774223 


-9.1182387633200 


-9.1182387832920 


(H) 


-5.3775999070685 


-5.3776042133885 
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-5.3775999070684 
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FIG. 3. Ground state of the dimensionless Hamiltonian ((34} with g = 15 (only positive values are shown). Left: u{q) = gPoo(g) 
[black solid line: wavefunction obtained directly in the momentum space - gray dashed line: wavefunction computed by 
Fourier transform of the wavefunction obtained directly in the configuration space]. Right: u{x) = a;7?.oo(s;) [black solid line: 
wavefunction obtained directly in the configuration space - gray dashed line: wavefunction computed by Fourier transform of 
the wavefunction obtained directly in the momentum space], h = 0.5 and N = 20 {h = 0.4 and A'^ = 20) for the computation 
in momentum (configuration) space. 



The wavefunctions produced by the LMM in both spaces have also been compared. A typical example is shown on 
Fig- 13 On the left part, the wavefunction obtained directly in the momentum space and the wavefunction computed by 
the Fourier transform of the wavefunction obtained directly in the configuration space are superposed. The situation 
is the opposite on the right part. In both cases, the agreement is very good for low values of the arguments, but the 
Fourier transform wavefunctions present large unphysical oscillation for larger values of the argument. The starting 
points of these oscillations can be rejected to high values of the argument in the asymptotic tail, but it is then 
necessary to work with greater values of N, typically N > 100 [Ij]. A good representation of the wavefunction in 



momentum space can be obtained with a small mesh only by working directly in this space. This an advantage of the 
LMM in momentum space. 

As in the nonrelativistic case, we will not try to study a "realistic" relativistic system. We will consider quite 
arbitrary, but convenient, values for the parameters. 



TOi = 7712 



1, 



1, 



(35) 



for which only one bound state exists (0 < E < 2m) . As we can see by examining Fig. 21 the behaviors of the eigenval- 
ues as a function of N and h is similar for nonrelativistic and semirelativistic Hamiltonians. But, the convergence is 
less good for semirelativistic systems, with shorter plateaus and larger variations around the plateaus. This situation 
is similar in configuration space [8|, |l3| . Various observables have been computed and are presented in Table [III A 
reasonable accuracy can be obtained even with a small mesh. 
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FIG. 4. Ground state eigenvalue E of ()22|l for the spinless Salpeter Hamiltonian with the Gaussian interaction considered (|35 
as a function of h and A'^. For A'^ = 10 and h = 1.0, the value of E is 1.8750 far above the range of the graph. 



TABLE II. Eigenvalues E and mean values of some observables for the ground state of the spinless Salpeter Hamiltonian 



with the Gaussian interaction U{r) 



-a exp I 



considered p5|) . as a function of A*' for h = 0.4. Computations in the 



momentum (Mom.) space are compared with accurate results obtained in the configuration (Conf.) space for A^ = 100 and 
h = 0.4. In order to check the mean values, (H) (which must be equal to E) is computed by 2{^/p^ + rn?) + {U{r)). 





Conf. 




Mom. 








Af = 10 


A- = 20 


A = 50 


E 


1.87098362 
1.3553804 


1.87044199 
1.3542724 


1.87100878 
1.3554650 


1.87098367 


Wp^+m^) 


1.3553807 


{Vp~^) 


3.991567 


3.981098 


3.992369 


3.991570 


(r) 


1.73375 


1.71171 


1.73551 


1.73376 


{U{r)} 


-0.8397772 


-0.8381094 


-0.8399212 


-0.8397777 


(H) 


1.87098362 


1.87043532 


1.87100880 


1.87098367 



B. Yukawa potential 

Some numerical tests are also performed using the Yukawa potential whose asymptotic behavior is very different in 
the configuration and momentum spaces [21| 



V{r) = -%xp(-6r) ^ VMk) = -^^jfp- 



(36) 



Using (dD), we can compute 
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Vi{p,p') 



npp' 



h"^ +p^ + p'' 



(37) 



where Qi{x) is a Legendre function of the second kind [35|- This type of potentials is used in many field of physics, 
even in hadronic physics. If the interaction in a quark-antiquark system or between two gluons in the vacuum can 
be simulated by a funnel potential (Coulomb+linear), the confinement vanishes inside a quark-gluon plasma above a 
critical temperature and the interaction could turn into a Yukawa potential [36[ . 
With a nonrelativistic kinematics, we can write 



52 



where e is the solution of the dimensionless Hamiltonian 



H = q'-g- 



with g — 



2/ia 



(38) 



(39) 



We choose to study the performance of the LMM in momentum space with the value g = 10. In this case, three 
bound states exist. A phenomenological approximate formula [37[ gives -E(„=o,i=o) — ~16.32, -E(i.o) = —0.65 and 



E, 



(0,1) 



-0.22. 



On Fig. \5\ we can see the behaviors of the ground state eigenvalue e oi H ([5^ . obtained by solving (|22p . as a 
function of N and h. The situation is at first sight similar to the case of the Gaussian potential, but the convergence 
is much slower for the Yukawa potential. For a mesh of 50 points, a plateau does not appear clearly, just a slowing 
down of the variation of e. For computations in the configuration space, we have checked that a small plateau is 
already present for iV = 20. Consequently, eigenvalues computed in momentum space for two different values of h will 
only coincide for large values of TV. About ten times more points are necessary for computations in momentum space 
to reach an accuracy similar to the one of computations in configuration space. It is illustrated in Table IIIII where 
eigenvalues and some observables are calculated with the LMM in both configuration and momentum spaces. Again, 
the mean value (H), which must be equal to e, is computed. The disagreement for the results in momentum space 
shows that the convergence is not so good than for the results in configuration space. We think that this peculiarity 
is linked to the asymptotic behavior of the Yukawa potential in the momentum space: Its decrease like a power-law 
(|36p is quite slow. The wavefunction in this space is then very extended and requires its computation at high values 
of p to be well described on its physical domain. 
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600 
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FIG. 5. Ground state eigenvalue e of (|22p for the dimensionless Hamiltonian (|39|) with g — 10. 



The wavefunctions produced by the LMM in both spaces have also been compared for the Yukawa potential on 
Fig. [S] With a small mesh of 20 points, the relative error on the ground state eigenvalues is around 10""'^° for 
the computation in configuration space, and around 10~^ for the computation in momentum space. Curiously, 
the wavefunction computed in momentum space is already very well reproduced and is even better than the Fourier 



11 



TABLE III. Eigenvalues e and mean values of some observables for H = q^ + U{x) with U{x) = ~ge ^ jx (|39p and g = 10. 
Computations in configuration (Conf.) and momentum (Mom.) spaces are performed with A*' = 200. The values of the 
parameter h are arbitrarily taken into a corresponding plateau. In order to check the mean values, {H) (which must be equal 
to e) is computed by (q'^) + {U{x)). 



n = 0, Z = 
Conf. Mom. 

{h = 0.02) {h = O.i 



n= 1, 1 = 

Conf. Mom. 

{h = 0.05) {h = 1.0) 



n = 0, 1 = 1 

Conf. Mom. 

{h = 0.05) {h = 0.5) 



e -16.340426 -16.340415 -0.6053933 -0.6053975 -0.205082327 -0.205082331 

{q^) 23.788977 23.788942 2.95238 2.95241 2.70792857 2.70792862 

{U{x)) -40.1294 -40.1200 -3.55778 -3.55743 -2.913010896 -2.913010877 

jH) -16.340426 -16.331047 -0.6053933 -0.6050217 -0.205082327 -0.205082257 




/ 



1.5 2.0 2.5 3.0 



FIG. 6. Ground state of the dimensionless Hamiltonian (|39|l with g = 10 (only positive values are shown). Left and right parts 
as in Fig. ^ h = 0.5, iV = 20 and e = -16.2066 {h = 0.05, A^ = 20 and e = -16.3404) for the computation in momentum 
(configuration) space. 

transform computed from the computation in configuration space with the same mesh. By increasing N, the beginning 
of the unphysical oscillations in the Fourier transforms can be rejected far in the asymptotic tail. 

We will briefly comment the Yukawa interaction with a semirelativistic kinematics. As in the nonrelativistic case, 
we did not try to study a "realistic" system. Quite arbitrary but convenient values are considered for the parameters, 



TOi = TO2 



16, 



1, 



(40) 



for which only one bound state exists {0 < E < 2m). Results about the convergence are presented in Fig. [T] Graphics 
are similar to the ones produced in the nonrelativistic case: Existence of plateaus for energy as a function of h] increase 
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FIG. 7. Ground state eigenvalue E of (|22[> for the spinless Salpeter Hamiltonian with the Yukawa interaction considered (140 |l . 
as a function of h and A. For A = 25 and h = 0.5, the value of E is 30.81 far above the range of the graph. 
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of the length of the plateau with N. But, the convergence is slower than in the nonrelativistic computations: N must 
be greater to achieve a quasi fiat plateau (Fig. [71 left) and to reach a convergence between eigenvalues computed with 
different values of h (Fig. [71 right). We have checked that it is also the case to obtain an agreement between observables 
computed in both spaces and to obtain better Fourier transforms. With a semirelativistic Hamiltonian, the kinetic 
energy increases as i/p^ not as p^ . So, higher values of the momentum can be reached by the wavefunction, and it 
can be more extended than with a nonrelativistic kinematics. This amplifies the problem mentioned in the previous 
section about the asymptotic behavior of the potential. 

V. CONCLUDING REMARKS 

The Lagrange-mesh method is a procedure to compute accurate eigenvalues and eigenfunctions of quantum equa- 
tions. Implemented at the origin in the configuration space, this technique requires only the computation of the 
potential at some mesh points Q, Q • This is due to the use of a Gauss quadrature rule with the fact that the basis 
functions satisfy the Lagrange conditions, that is to say they vanish at all mesh points except one. 

Using this very special property, we have shown that the Lagrange-mesh method can be adapted to be used directly 
in momentum space with a nonlocal potential (but local in the configuration space). In this case, nonrelativistic 
and semirelativistic kinetic parts are treated with the same manner. Moreover, it is possible to treat interactions 
which present discontinuities in configuration space. The wavefunction in momentum space is directly obtained under 
the form of expansion coefficients for the Lagrange functions. Using only these coefficients, it is easy to compute 
the wavefunction in configuration space by Fourier transform, as well as mean values of p-dependent or r-dependent 
operators. Though the technique is not variational, the convergence can be checked: Eigenvalues present a plateau as 
a function of the sole nonlinear parameter of the method; eigenvalues computed with different nonlinear parameters 
in a plateau tends towards the same value as the number of mesh points increases. 

The method has been tested with two potentials. With a Gaussian interaction and a nonrelativistic kinematics, 
a good accuracy can be obtained for eigenvalues and observables with a very small mesh, as for the method in 
configuration space: 20 points are probably sufficient for most of the physical applications. For a Yukawa interaction, 
it is necessary to use more points than for the Lagrange-mesh method in the configuration space, typically 10 times 
more. For both potentials: i) The convergence is slowed down for a semirelativistic kinematics, as for the method in 
configuration space. But a good accuracy can nevertheless be reached [B|; H) It is necessary to use a large mesh to 
obtain a correct Fourier transform of the wavefunctions, as for the computations in configuration space [13|- Hi) A 
good wavefunction in momentum space can obtained with a small mesh by a direct computation in momentum space. 
It does not present the unphysical oscillations of the Fourier transform of the wavefunctions in configuration space, 
which can only be eliminated by using a large mesh. 

The purpose of this paper is to test the feasibility and the reliability of the Lagrange-mesh method in momentum 
space. We think that it can be safely applied to physical problems, since convergence tests exists and a good accuracy 
can always be obtained. Even if several hundreds of points are necessary to treat some potentials, eigenvalues can be 
rapidly computed. The method is very easy to implement once the interaction is known in momentum space, and it 
can be particularly useful if the kinetic part T{p^) is an unusual function of the momentum, since its corresponding 
matrix is diagonal and easy to compute. This kind of situation appears in hadronic physics where quarks or gluons 
can be considered with a momentum dependent mass which can be very complicated to define in the configuration 
space [28t - t3C| , e.g. T(p^) = ^/jp^+rn?(jp) . For these systems, the dominant interaction can be simulated by a linear 
potential in the configuration space. This type of potential is a priori highly singular in the momentum space, but it 
exists in a distributional sense and can be used in this context ISa . 
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